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Abstract 

Model uncertainties or simulation uncertainties occur in mathematical modeling of 
multiscale complex systems, since some mechanisms or scales are not represented (i.e., 
"unresolved") due to lack in our understanding of these mechanisms or limitations in 
computational power. The impact of these unresolved scales on the resolved scales 
needs to be parameterized or taken into account. A stochastic scheme is devised to 
take the effects of unresolved scales into account, in the context of solving nonlinear 
partial differential equations. An example is presented to demonstrate this strategy. 
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1 Introduction 



Mathematical models for scientific and engineering systems often involve with some un- 
certainties. We may roughly classify such uncertainties into two kinds. The first kind of 
uncertainties may be called model uncertainty. They involve with physical processes that 
are less known, not yet well understood, not well-observed or measured, and thus difficult 
to be represented in the mathematical models. 

The second kind of uncertainties may be called simulation uncertainty. This arises 
in numerical simulations of multiscale systems that display a wide range of spatial and 
temporal scales, with no clear scale separation. Due to the limitations of computer power, 
at present and for the conceivable future, not all scales of variability can be explicitly 
simulated or resolved. Although these unresolved scales may be very small or very fast, 
their long time impact on the resolved simulation may be delicate (i.e., may be negligible 
or may have significant effects, or in other words, uncertain). Thus, to take the effects of 
unresolved scales on the resolved scales into account, representations or parameterizations 
of these effects are desirable. 

These uncertainties are sometimes also called unresolved scales, as they are not repre- 
sented or not resolved in modeling or simulation. Model uncertainties have been considered 
in, for example, [lQj [HI CEl [2J EU [T71 [26j ETJ [28] and references therein. Works relevant for 
parameterizing unresolved scales include [15j 021 EH S 001 EH [331 S] > among others. 

In this paper we consider an issue of approximating model uncertainty or simula- 
tion uncertainty (unresolved scales) by stochastic processes, and then devise a stochastic 
scheme for such approximations. We first recall some basic facts about fractional Brown- 
ian motion (fBM) in Sj2j Then we discuss model uncertainty and simulation uncertainty 
in Sj3]and £JH respectively. Finally, we present an example in f|5] demonstrating our result. 
This example involves approximating subgrid scales via correlated noises, in the context 
of large eddy simulations of a partial differential equation. 

2 Fractional Brownian motion and colored noise 

We discuss a model of colored noise in terms of fractional Brownian motion (fBM), includ- 
ing a special case which is white noise in terms of usual Brownian motion. The fractional 
Brownian motion B H (t), indexed by a so called Hurst parameter H G (0,1), is a gen- 
eralization of the more well-known process of the usual Brownian motion B(t). It is a 
centered Gaussian process with stationary increments. However, the increments of the 
fractional Brownian motion are not independent, except in the usual Brownian motion 
case (H = \). For more details, see [25| I2B1 [HI 120} [38] . 

Definition of fractional Brownian motion: For H G (0, 1), a Gaussian process B H (t), 
or fBM it), is a fractional Brownian motion if it starts at zero B H (0) = 0, a.s., has mean 
zero E[B H (t)} = 0, and has covariance E[B H (t)B H (s)} = ±(\t\ 2H + \s\ 2H - \t - s\ 2H ) for 
all t and s. The standard Brownian motion is a fractional Brownian motion with Hurst 
parameter B =\. 
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Some properties of fractional Brownian motion: A fractional Brownian motion B (t) 
has the following properties: 

(i) It has stationary increments; 

(ii) When H = 1/2, it has independent increments; 

(hi) When H ^ 1/2, it is neither Markovian, nor a semimartingale. 

We use the Weierstrass-Mandelbrot function to approximate the fractional Brownian 
motion. The basic idea is to simulate fractional Brownian motion by randomizing a 
representation due to Weierstrass. Given the Hurst parameter H with < H < 1, we 
define the function w(t) to approximate the fractional Brownian motion: 

oo 

w(ti) = Cji J " >\u{2-r J i, ■ (I,) 

j=-oo 

where r = 0.9 is a constant, Cj's are normally distributed random variables with mean 
and standard deviation 1, and the dj's are uniformly distributed random variables in the 
interval < dj < 2ir. The underlying theoretical foundation for this approximation can 
be found in [31, 22] , Figures 1 and 2 show a sample path of the usual Brownian motion 
(i.e., H = ^), and fractional Brownian motion with Hurst parameter H = |, respectively. 



2 




-2 1 1 1 1 1 1 

5 10 15 20 25 30 

t 

Figure 1: A sample path of Brownian motion B{t) 

3 Model uncertainty 

We consider a spatially extended system modeled by a partial differential equation (PDE): 

u t = Au + N(u), (1) 
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Figure 2: A sample path of fractional Brownian motion B 



(t), with H = 0.75 



where A is a linear (unbounded) differential operator, and iV is a nonlinear function of 
u(x,t) with x € D and t > 0, and satisfies a local Lipschitz condition. In fact, N may also 
depend on the gradient of u. 

If this (deterministic) model is accurate, i.e., its prediction on the field u matches 
with the observational data u on a certain period of time [0,T], then there is no need for 
a stochastic approach. However, when the prediction u deviates from the observational 
data u, we then need to modify the model (JTJ) . In this case, the observational data u may 
be thought to satisfy a modified model: 



where the model uncertainty F{u) is usually a fluctuating (i.e., random) process, as the 
observational data u is so (i.e., has various samples or realizations). 

The model discrepancy or model uncertainty F(u) may have various causes, such as 
missing physical mechanisms (not represented in the deterministic model (HJ). Sometimes, 
the model uncertainty F(u) is smaller in magnitude than other terms in the model (J2|) 
and thus is often ignored in the deterministic modeling. However, being small and being 
fluctuating may not necessarily imply that its impact on the overall system evolution to 
be small pQ. To take this impact into account, we would like to model or approximate 
F(u) by a stochastic process. 

We first calculate the model uncertainty F(u) via observational data u. By discretizing 
(|2|) and using data samples for u, we obtain (discretized) samples for F. 

The time correlation may then be calculated using the samples of F. If the time 
correlation scale is significantly shorter than the time scale for the field u, we may ignore 




(2) 



the time correlation and thus approximate F by the following stochastic process containing 
a (uncorrelated) white noise, for example: 

F = f + auB t , (3) 

where / = EF is the mean of F (computed from data), Bt is the usual Brownian motion 
(reviewed in ^5] below) and the deterministic noise intensity a may depend on space. Here 
a may be computed via stochastic calculus, especially the Ito isometry, as follows. 

F -EF = a uB t , 

r T 

EF)dt] 2 = o 2 E[ udB t ] 2 , 
Jo 

r T 

EF)dt} 2 = a 2 E u 2 (x,t)dt, 
Jo 



E[f (F - 
Jo 

E[f (F - 
Jo 



Thus we obtain 



a 



\ 



E[£(F-EF)dt}< 
E J^u 2 (x,t)dt 



(4) 



With this approximation, we obtain the following stochastic partial differential equation 
(SPDE) as a modified model for the original deterministic model ([!]):[/ ~ u 



U t = AU + N(U) + f(U) + a uB t . 



(5) 



In general, the model uncertainty F may be better approximated by correlated noise 
via fractional Brownian motion. Since the procedure is similar, we will demonstrate this 
in the next section when we discuss simulation uncertainty. 



4 Simulation uncertainty 

This section deals with simulation uncertainty, i.e., stochastically parameterizing the ef- 
fects of the unresolved scales on the resolved scales. We consider this issue in the context 
of large eddy simulations (LES) of a nonlinear partial differential equation with memory. 

In large eddy simulations of fluid or geophysical fluid flows [331 0], the unresolved 
scales appear as the so-called subgrid scales (SGS). The SGS term appears to be highly 
fluctuating ("random"); see the Figure 1 in [23]. Partially motivated by this, stochastic 
parameterizations of subgrid scales have been investigated in fluid, geophysical and climate 
simulations, based on physical or intuitive or empirical arguments. Another, perhaps more 
important, motivation for applying stochastic parameterizations of subgrid scales is to 
induce the desired backward energy flux ("stochastic backscatter" ) in fluid simulations 

PI EH ESI- 

We present one stochastic parameterization scheme of the subgrid scale term in the 
large eddy simulation of a nonlinear partial differential equation with an extra memory 
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term, which is in fact a nonlinear integro-partial differential equation. The approximation 
scheme is based on stochastic calculus involving with a fractional Brownian Motion, and 
the "parameter' to be calculated is a spatial function, which is derived using Ito stochastic 
calculus. 



u t = Au + F(u), (6) 

where A is a linear differential operator, and F is a nonlinear function of u(x, t) with 
x € .Dand t > 0, and satisfies a local Lipschitz condition. We investigate stochastic 
parameterizations of unresolved scales in the context of large eddy simulations of the 
above system. 

The idea of large eddy simulation is to split the flow into a local, spatial mean (or 
average) and a fluctuation about the that mean. The mean u is defined by filtering or 
mollification (convolution with an approximate identity). The goal is to predict the mean 
accurately. This is widely believed possible based on the idea that since fluctuations have 
random character, their average effects on the mean notion can successfully be medelled. 

To filter the solution, we pick a filter. Many different ones are commonly used. To fix 

the ideas, in this paper, we use Gaussian filter as in [4j, Gg(x) = -^%zz * T , where 5 > is 
the filter size and the filter is such that: (i) u* Gs is infinitely differentiable in space and, 
(ii) u * Gs — > u as 5 — > in L 2 (D). Here and hereafter u * Gg = f D u(y,t)Gg(x — y)dy or 
the over bar u denotes convolution. 

Remark 1. The mean u is a weighted average of u about the point x. As 5 — > 0, the 
points near x are weighted more and more heavily, so u * Gs — >• u as 5^0 in L 2 {D). 

Using the fact that convolution commutes with differentiatian, we get the space-filtered 
system: 



u t = Au + F(u), 

or 

u t = Au + F(u) + R(u,u), (7) 



where the subgrid scale term R(u,u) := F(u) — F(u). Since generally F{u) ^ F{u), the 
usual parameterization or closure problem of the large eddy simulation has arisen. Due 
to inaccurate (uncertain) initial conditions or boundary conditions, R(u, u) is a correlated 
fluctuating process [MIE]; depending on samples wina suitable sample space Q. We thus 
would like to approximate the subgrid scale term R(u, u) by a stochastic process with a 
correlated (i.e., colored) noise component, for example: 

dB H 

R = f(u) + a(x)^ r , (8) 

dB H 

where — -jt- is a colored noise (generalized time derivative of a fractional Brownian motion; 
reviewed in f|5] below), and 

f(u) = ER, (9) 



is the mean component of the subgrid scale term R. Moreover, the noise intensity cr(x) is 
a non-negative deterministic function to be determined from fluctuating SGS data R. The 
subgrid scale term R may be inferred from observational data (see |29l 150] for relevant 
information for subgrid scales in Navier-Stokes equations), or from fine mesh simulations. 

Note that a is to be calculated or estimated from the fluctuating SGS data R, either 
from observation or from fine mesh simulations. So this is an inverse problem. As in 
usual inverse problems [36], the stochastic parameterizations for the SGS term R is not 
unique. What we proposed above is merely an example. This offers an opportunity for 
trying various stochastic parameterization schemes, much as one uses various smoother 
functions (e.g., polynomials or Fourier series) to approximate less regular functions or data 
in deterministic approximation theory. 

To estimate the unknown parameter (function) <j(x), we start with the following rela- 
tion: 

R-ER = a{x)^j-. (10) 
Taking time integral over a computational interval [0, T] on both sides, we obtain 

r-T ,-T 



[ [R- ER]dt = I (r(x)dBt 1 = a(x)B$ . 
Jo Jo 



Therefore, taking mean-square on both sides 

r-T 



E( / [R-ER]dtf = a 2 {x)T 2 
Jo 

Thus an estimator for cr(x) is 



a{x) = ^JE( J T [R-ER]dty , (11) 



which can be computed numerically. 

By the stochastic parameterization ([5} on the SGS term R, with / determined from 
and a from (jlip . the LES model ([7]) becomes a stochastic partial differential equation 

(SPDE) for the large eddy solution U ~ u: 

dB H 

U t = AU + N(U) + f(U) + *{x)-£-, (12) 
with the appropriately filtered boundary condition and filtered initial condition. 



5 An example 

We present a specific example of stochastic modeling of simulation uncertainty of subgrid 
scales, in the context of large eddy simulations. We consider the following nonlinear partial 
differential equation with a memory term (time-integral term) [6]: 

/"* 1 

ut = u xx + u-u 3 + — — n ,u(x,s)ds, (13) 

Jo 1 + 1*- sf 
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under appropriate initial condition u(x,0) = uq(x) and boundary conditions u{ — l,t) = 
a, u(l,t) = b with a, b constants, on a bounded domain D : —1 < x < 1. Here (3 is 
a positive constant. This model arises in mathematical modeling in ecology [52], heat 
conduction in certain materials [111 [T9] and materials science [9] [19] . The time-integral 
term here represents a memory effect depending on the past history of the system state, 
and this memory effect decays polynomially fast in time. 

The large eddy solution u is the true solution u looked through a filter: i.e., through 
convolution with a spatial filter G$(x), with spatial scale (or filter size or cut-off size) 
5 > 0: 

u(x,t) := u* G s = / u(y,t)G s (x - y)dy. 
Jd 



,2 



In this paper, we use a Gaussian filter as in [5], G&(x) = ^p-e s 3 ". 
On convolving (fl~3|) with Gs, the large eddy solution u is to satisfy 



1 



u t = u xx + u - u 3 + I 1 + > t _ S \/3 U ( X ' s ) ds > 







or 



, ft 1 

u t = u xx + u-u + I — — ^ u(x,s)ds + R(x,t), (14) 

Jo l + l*- »r 

where the remainder term, i.e., the subgrid scale (SGS) term R(x,t) is defined as 

R(x,t) := (u) 3 - <y). (15) 

We can write u = u + u' with u the large eddy term and u' the fluctuating term. Note 
that u = u — v! . So the SGS term R(x,t) involves nonlinear interactions of fluctuations 
u' and the large eddy flows. Thus R(x,t) may be regarded as a function of u and u': 
R := R(u,u'). 

The leads to a possibility of approximating R(x, t) by a suitable stochastic process 
defined on a probability space (0,JF, P), with u S fi, the sample space, a— field .F and 
probability measure P. This means that we treat R data as random data as in [25] , which 
take different realizations, e.g., due to fluctuating observations or due to numerical simu- 
lation with initial and boundary conditions with small fluctuations. In fluid or geophysical 
fluid simulations, the SGS term may be highly fluctuating and time-correlated [24] . and 
this term may be inferred from observational data [29] 130]. or from fine mesh simulations. 

This further suggests for parameterizing the subgrid scale term R(x,t) as a time- 
correlated or colored noisy term. The increments of fractional Brownian motion are corre- 
lated in time and hence its generalized time derivative Bp is used as a model for colored 
noise. In the special case H = ^, we have the white noise Bf. Thus we parameterize the 
subgrid scale term R(x,t), which is time-correlated, by colored noise Bp as follows: 

R( x ,t)=f(u) + a(x)^f, (16) 
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where 



f(u) =MR(x,t), (17) 

is the mean component of the subgrid scale term R(x,t). Moreover, the noise intensity 
&{x) is a non-negative deterministic function to be determined from fluctuating SGS data 
R. The subgrid scale term R(x, t) may be inferred from observational data [291 130] . or from 
fine mesh simulations as we do here. We represent the mean component f(u) in terms of 
the large eddy solution u. The specific form for / depends on the nature of the mean of 
R. Here we take f(u) = ao + a\u + a2U 2 + a^u 3 , where coefficients aj's are determined via 
data fitting by minimizing J^Jao + aiu + a2U 2 + o,^u 3 — Ei?(x, t)] 2 dxdt. Moreover, we 
take Bf as a scalar fractional Brownian motion. 

Note that a is to be calculated or estimated from the fluctuating SGS data R, either 
from observation or (in this paper) from fine mesh simulations; see detailed discussions in 
[241 \7\. So this is an inverse problem. As in usual inverse problems [36], the stochastic 
parameterizations for the SGS term R is not unique. This offers an opportunity for 
trying various stochastic parameterization schemes, much as one uses various smoother 
functions (e.g., polynomials or Fourier series) to approximate less regular functions or data 
in deterministic approximation theory. 

To estimate the unknown parameter (function) cr(x), we start with (|161) -(fT7]) to get 
the following relation: 

dB H 

R(x,t) -ER(x,t) = a(x)—^-. (18) 
Taking time integral over a computational interval [0, T] on both sides, we obtain 

r T r T 

/ [R(x, t) - ER(x, t)]dt = / a(x)dBf I = a{x)B% . 
Jo Jo 

Therefore, taking mean-square on both sides, 

E(/ [R(x,t) -ER(x,t)]dt) 2 = a 2 (x)T 2H . 
Jo 

Thus an estimator for cr(x) is 



a(x) = ^ w JE(j\R(x,t)-ER(x,t)}dt) 2 , (19) 



which can be computed numerically. 

By the stochastic parameterization (j!6j) on the SGS term R, with / determined from 
(|17p and a from (|19p . the LES model (|14p becomes a stochastic partial differential equation 
(SPDE) for the large eddy solution U ~ u: 

ft i dB H 

U t = U xx + U -U 3 + / — — ^ U(x, s)ds + f(U) + cj(x)—^, (20) 

J Q 1 + \t - s\P dt 
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with boundary conditions U(—l,t) = a, U(l,t) = b and filtered initial condition 



U(x, 0) = uq(x). 



(21) 



Numerical Experiments: 

We use a spectral method to solve nonlinear system (|13p and (|20p numerically. For 
more details, please see [37]. We take the following initial and boundary conditions: 



Fine mesh simulations of the original system with memory (|13p are conducted to gener- 
ate benchmark solutions or solution realizations, with initial conditions slightly perturbed; 
see Fig. 3. These fine mesh solutions u are used to generate the SGS term R defined in 
(1151) at each time and space step. The filter size used in calculating R is taken as 5 = 0.01. 
The mean / is calculated from (|17p via cubic polynomial data fitting (as discussed in the 
last section), and parameter function a(x) is calculated as in ([T9]h The stochastic LES 
model ()20p is solved by the same numerical code but on a coarser mesh. Note that a four 
times coarser mesh simulation with no stochastic parameterization for the original system 
(fT3l) does not generate satisfactory results; see Fig. 4. The stochastic LES model (|20P is 
then solved in the mesh four times coarser than the fine mesh used to solve the original 
equation (I13p , The stochastic parameterization leads to better resolution of the solution 
as shown in Fig. 5. As in [7J, it can be shown that when two stochastic parameterization 
terms are close in mean-square norm on finite time intervals, the solutions are also close 
in the same norm. 



u(x,0) = uq = 0.53x — 0.47sm(1.57rx), u(— 1, i) 



1, «(l,t) = l 



=. o- 
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100 




20 
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Figure 3: Solution to the original system on a fine mesh, ut = u xx + u — W 
= 2, mesh size Ax = 0.001. 




o l+|f-« 



■p u(x, s)ds, 
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Figure 4: Solution to the original system with NO stochastic parametrization on the mesh four 
times coarser than the mesh used in Fig. 3, u t = u xx + u — u 3 + J* 1+ u_ s w u(x, s)ds, (3 = 2. 




Figure 5: Solution to LES model with stochastic parametrization on the mesh four times coarser 
than the mesh used in Fig. 3, U t = U XX + U - U 3 + f* „ U(x, s)ds + f(U) + a{x)B^, (3 = 2. 
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